function [ A, S, B, out] = SoftImputeALS( D, lambda, maxR, para )

maxIter = para.maxIter;
decay = para.decay;
tol = para.tol;

m = size(D, 1);
n = size(D, 2);

[row, col, data] = find(D);

R = randn(n, 1);
U = powerMethod( D, R, 5, 1e-6);
[~, S, V] = svd(U'*D, 'econ');
lambdaMax = S;

A = eye(m, maxR);
A(:, 1) = U;
B = eye(n, maxR);
B(:, 1) = V;

clear U V R;

spa = D;
obj = zeros(maxIter, 1);
Time = zeros(maxIter, 1);
RMSE = zeros(maxIter, 1);
t = tic;
for i = 1:maxIter
    lambdai = abs(lambdaMax - lambda)*(decay^i) + lambda;
    
    temp = partXY(A', B', row, col, length(data));
    temp = data - temp';
    spa = setSval(spa, temp, length(temp));
    
    Ai = B'*B;
    Ai = Ai + lambdai*eye(size(Ai));
    Ai = pinv(Ai);
    A = spa*(B*Ai) + A*(B'*B)*Ai;
    
    temp = partXY(A', B', row, col, length(data));
    temp = data - temp';
    spa = setSval(spa, temp, length(temp));
    
    Bi = A'*A;
    Bi = Bi + lambdai*eye(size(Bi));
    Bi = pinv(Bi);
    B = spa'*A*Bi + B*(A'*A)*Bi;
    
    obji = sum(temp.^2)/2;
    obji = obji + (lambda/2)*sum(A(:).^2);
    obji = obji + (lambda/2)*sum(B(:).^2);
    obj(i) = obji;
    
    if(i == 1)
        delta = inf;
    else
        delta = abs(obj(i - 1) - obj(i));
    end
    
    % testing performance
    Time(i) = toc(t);
    fprintf('iter %d; obj:%.3d (%.3d) \n', i, obji, delta);
    if(isfield(para, 'test'))
        tempS = eye(size(A, 2), size(B, 2));
        RMSE(i) = MatCompRMSE(B, A, tempS, ...
            para.test.row, para.test.col, para.test.data);
        fprintf('RMSE %.2d \n', RMSE(i));
    end
    
    if(delta < tol)
        break;
    end
end

[A, RA] = qr(A, 0);
[B, RB] = qr(B, 0);

[U, S, V] = svd(RA*RB', 'econ');
S = diag(S);
nnzS = sum(S > 1e-5);
S = S(1:nnzS);
S = diag(S);

out.Rank = nnz(S);

A = A*U(:, 1:nnzS);
B = B*V(:, 1:nnzS);

out.obj = obj(1:i);
out.RMSE = RMSE(1:i);
out.Time = Time(1:i);

end

